The two-dimensional random-bond Ising model, free fermions and the network model 
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We develop a recently-proposed mapping of the two-dimensional Ising model with random exchange 
(RBIM), via the transfer matrix, to a network model for a disordered system of non-interacting 
fermions. The RBIM transforms in this way to a localisation problem belonging to one of a set of 
non-standard symmetry classes, known as class D; the transition between paramagnet and ferromag- 
net is equivalent to a delocalisation transition between an insulator and a quantum Hall conductor. 
We establish the mapping as an exact and efficient tool for numerical analysis: using it, the compu- 
tational effort required to study a system of width M is proportional to M 3 , and not exponential 
in M as with conventional algorithms. We show how the approach may be used to calculate for the 
RBIM: the free energy; typical correlation lengths in quasi-one dimension for both the spin and the 
disorder operators; even powers of spin-spin correlation functions and their disorder-averages. We 
examine in detail the square-lattice, nearest-neighbour ± J RBIM, in which bonds are independently 
antiferromagnetic with probability p, and ferromagnetic with probability 1 — p. Studying tempera- 
tures T > 0.4J, we obtain precise coordinates in the p — T plane for points on the phase boundary 
between ferromagnet and paramagnet, and for the multicritical (Nishimori) point. We demonstrate 
scaling flow towards the pure Ising fixed point at small p, and determine critical exponents at the 
multicritical point. 



I. INTRODUCTION 

The two-dimensional Ising modelBi has been a basic 
prototype in the theory of phase transitions for over half 
a century. A central factor in its importance has been 
its equivalence to a system of non-interacting fermions, 
as set out by Schultz, Mattis and LiebH in their well- 
known reformulation of Onsager's solution. The two- 
dimensional Ising model has naturally also been a test- 
bed for studies of the effect of quenched disorder on phase 
transitions, and the equivalence between the spin system 
and free fermions continues to hold in the presence of ran- 
domness in exchange interactions. Iruthis paper we build 
on recent work bv-Cho and Fisherfm and by Gruzberg, 
Read and LudwigtlQ to establish the correspondence in a 
form suitable for numerical analysis, and use it to study 
the square-lattice, random-bond Ising model (RBIM). 

The consequences for the two-dimensional Ising model 
of weak randomness in exchange interactions are rather 
well understood, following analytical calculations based 
on the fermionic formulation by Dotsenko and DotsenkoO 
and others :ErO weak disorder is marginally irrelevant 
in the renormalisation group sense, and the thermally- 
driven transition from the paramagnet to the ferromag- 
net survives with only logarithmic modifications to the 
critical behaviour of the pure system. By contrast, strong 
disorder has more dramatic effects. A convenient choice 
is to consider exchange interactions with fixed magni- 
tude which are independently ferromagnetic or antiferro- 
magnetic, with probabilities 1 — p and p respectively Lu. 
this case, it is known from a variety of approachesc3~E2l 
that the Curie temperature is depressed with increasing 
p, reaching zero at a critical disorder strength, p c . More- 
over, while the scaling flow at the transition is controlled 



for small p by the critical fixed point of the pure sys- 
tem, at larger p it is determined by a disorder-domiaatcd 
multicritical point, known as the Nishimori point .llailij 

Most numerical studies of— tbe RBIM have used ei- 
ther Monte Carlo simulatianaiiEj or transfer matrix cal- 
culations in a spin basis.EIrEj Fermionic formulations 
of the Ising model nevertheless have two great poten- 
tial advantages: they can avoid the statistical sam- 
pling errors of Monte Carlo simulations; and also, if 
implemented using the transfer matrix, they can avoid 
the exponential growth in transfer matrix dimension 
with system width that occurs if this matrix is writ- 
ten in a spin basis. Pioneering steps in the fkst of 
these directions have been rtakfin by Blackmar£j and 
collaborators]^ and others,E3cj using the solution of 
the two-dimensional Ising model via a PfaffianH to ex- 
press statistical-mechanical quantities in terms of spec- 
tral properties of the associated matrix. Their work 
makes a link between the RBIM and localisation prob- 
lems, since the matrix allied to the Pfaffian is essentially 
a tight-binding Hamiltonian on the lattice of the underly- 
ing Ising model, with random hopping arising from ran- 
dom exchange. An alternative route from the RBIM to 
a localisation problem has been proposed by Cho and 
Fisher hO starting from two copies of the transfer matrix 
for an Ising model, each expressed in terms of Majorana 
fermions and combined to form Dirac fermions, they ar- 
rive at a version of the network model similar to that 
introduced as a description for the integer quantum Hall 
plateau transition]^ though with a distinct symmetry. 

Viewed localisation problem, the paramagnetic 

and ferromagnetic phases of the RBIM translate to two 
insulating phases with Hall conductance differing by one 
quantum unit, while the Curie transition maps to a 
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version of the quantum Hall plateau transition. This 
transition, and indeed the insulating phases, belong to 
a non-standard symmetry class for localisation, classi- 
fied in work by Altland and Zirnbauero and known as 
class D. The match between behaviour expected in the 
RBIM and that anticipated for two-dimensional localisa- 
tion proMi^s-in-class D has been the subject of recent 
discussion.QiJ'E3 cil A particular difficulty has been to rec- 
oncile the fact that, generically, a third, metallic phase 
is possible in the localisation problem, in addition to the 
two insulating phases, while the RBIM in two dimensions 
is expected to display only two phases. The resolution 
which has emergedtrcil is that symmetry alone is not suf- 
ficient to determine the phases that appear, and that in 
the specific disordered conductor equivalent to the RBIM 
no metallic phase arises. 

The work we describe here builds on Cho and Fisher's 
ideas, which must be extended in several ways to provide 
a precise and practical treatment of the RBIM. First, the 
approach described in Rcf. [| proceeds from the RBIM 
via a continuum limit, which is rediscretised to obtain a 
network model. In order to find an explicit relationship 
between parameters in the two systems, it is necessary 
instead to carry out the mapping directly on a lattice 
model. Doing so, as described by Cho in her thesisa and 
by Gruzberg, Read and Ludwig in Refs. ^ and |?], one ar- 
rives at a network model different in detail to. that stud- 
ied in Ref. ^, and with different behaviour.cd Second, a 
proper treatment of the RBIM in cylindrical geometry 
requires an appropriate choice of boundary conditions in 
the network model, which has not previously been con- 
sidered. Third, to calculate thermodynamic quantities, 
typical correlation lengths, spin and disorder correlation 
functions for the RBIM using the network model formula- 
tion, it is necessary to map from fermions back to spins, 
as outlined in Refs. || and [?] and as we describe here. 
A feature of interest which emerges from our analysis is 
a topological distinction between the paramagnetic and 
ferromagnetic phases as represented in terms of fermions, 
similar to that discussed recently for other systems from 
symmetry class D.QI13 Finally, an important technical as- 
pect of the work we present here is that numerical trans- 
fer matrix calculations for localisation problems in the 
symmetry class we are concerned with require for numer- 
ical stability a modification of the standard algorithm, as 
first discussed in Ref |4l]. 

As a numerical approach to the RBIM, the method we 
describe has two main limitations. One arises because 
the Dirac fermions of the network model are built from 
two copies of an Ising model. As a result, it turns out to 
be straightforward to calculate even powers of spin cor- 
relations functions, and their disorder-averages, but not 
practical to calculate odd powers. The other stems from 
the fact that Boltzmann factors which enter the network 
model become large at low temperatures, making the zero 
temperature limit inaccessible. 

The rema inder of the paper is organised as follows: 
In Sec. II A and Sec. II B we outline the Jordan- Wigner 



fermionisation of the spin transfe r ma trix and the map- 
ping to a network model. In Sec. II C we discuss bound- 



ary conditions across the system in network model lan- 
guage and the subsequent rules for constructing the spin 
transfer matrix from the fermion transfer matrix. In 
Sec. HI and Appendix^ we review the numerical algo- 
rithm that we employ in the network model transfer ma- 
trix calculations and set out how statistical-mechanical 
quantities are obtained from the fermion description. In 
Sec. IV we present numerical results on the ± J RBIM. 
The system sizes we study (transverse width M = 8 — 256 
spins) are significantly larger than what was previously 
possible. We focus on critical behaviour at the Nishi- 
mori point, for which we determine the coordinate p c = 
0.1093±0.0002. We calculate the critical exponents v and 
vt, describing the divergence of the correlation length as 
the Nishimori point is approached along the Nishimori 
line and the phase boundary, respectively. Using large 
system sizes we find v ^._L50 ± 0.03, in disagreement 
with previous estimates,t§E3 and, with wider confidence 
limits, v T = 4.0 ± 0.5. 



II. TRANSFER MATRIX 

A. Ising model transfer matrix 

We consider the nearest neighbour Ising model on a 
square lattice in two dimensions. The partition function 
Z for a such a,system on a strip of length L and width M 
can be written!! in terms of a product of transfer matrices. 
Introducing integer coordinates, n and i, as illustrated in 
Fig.0, one has 



Z = A Tr 



T X T 2 ■ ■ ■ T n ■ ■ ■ T L 



(2.1) 



with T n = V n x H n and 



H n = cxp (- Y^Li [<,i< ] 
V n = exp i K n,icrfcrf +1 ] 



(2.2) 



where the cr's are Pauli matrices and 



A 



iln [tanh (3Jh(n,i)], 



ULiUZi V2[sinh2< 



(2.3) 



Here, is the reduced coupling strength at inverse tem- 
perature (3 between the ith and (i + l)th spin in the ver- 
tical direction of Fig.|] on the nth slice, and k* ■ is the 
Kramers- Wannier dual value of the corresponding bond 
strength in the horizontal direction. For the rest of the 
paper the labels v and h on the bond strengths are re- 
dundant, since all horizontal bond strengths (and only 
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those) appear as dual valu es, i dentified with an asterix. 
We take o\ l+1 = of in Eq. ( |2.2| ) so that boundary condi- 
tions across the strip are controlled by the set of interac- 
tions strengths K n ,M- For convenience we introduce the 
notation T(k, I) = IlL=fc and f° r brevity we use T to 
denote either T(k, I) or T n . 



M 



i+1 



i+2 



J v (n,i-1) 



i(n,i) 



J v (n,i) 



J(n,i + 1) 



n+1 



n+2 



FIG. 1. For the square lattice Ising model we adopt the 
convention that a pair (n, i) labels one spin with two asso- 
ciated bonds, one horizontal (to the right) and one vertical 
(downwards) . 



Following Schultz, Mattis and LiebEl the operators 
H n and V n can be written, using the Jordan- Wigner- 
transformation, as functions of fermionic operators. In- 
troducing the fermion annihilation and creation opera- 
tors C{ and Cj, the spin operators become 



of =eK P (i7r£5=ictC i ) (Cj+d), 
af = 2dd - 1 . 



(2.4) 



After Jordan- Wigner transformation, H n and V n read 
H n = exp(-2£^ 1 <AC\C^\]) , 

V n = cxp ( Efir 1 «n,i [C\ ~ Ci] [Cj +1 + C i+1 ] (2.5) 
-Kn,Me™ Na [Cl - C M ][C\ + d] 



with Nc = ^2fLi C[Ci, the number operator. A famil- 
iar feature of the transfer matrix in fermionic language 
is that it does not conserve Nq, since V n includes terms 
which create and annihilate fermions in pairs. Such a 
structure is reminiscent of Bogoliubov-de Gennes Hamil- 
tonians arising in the mean-field description of super- 
conductors. It has the consequence that, to diagonalise 
the transfer matrix for a translationally invariant Ising 



model, one uses Fourier transformation followed by Bo- 
goliubov transformation. For the RBIM without trans- 
lational invariance, the transformation that diagonalises 
the transfer matrix is disorder-dependent, and one must 
follow a different route to make progress. 

In place of diagonalisation, the objective for the RBIM 
is to write the transfer matrix in terms of Dirac fermions 
whose number is conserved under., its action. The nec- 
essary steps are well-establishecOEl and have been set 
out in the present context by Cho and Fisher ,□ Cho,El 
and Gruzberg, Read and LudwigQ. First, because of the 
form of Eq. ( [2.5| ), it is natural to decompose the complex 
(Dirac) fermions into real and imaginary parts, introduc- 
ing real (Majorana) fermions £c and rjc- Suppressing the 
site index one can write 



C f = -^=[& + iVc] 



(2.6) 



where £c and f]c anticommute and satisfy = £c, 

Vc = and {£ci,£cy} = {vct,Vcj} = k r Next, in or- 
der to return to Dirac fermions, one introduces a second, 
identical copy of the Ising model. We represent the sec- 
ond copy using the Dirac fermions D and D' , in analogy 
to the C and C\ and employ the Majorana decomposi- 
tion D=[£ D - ir] D ]/V2 and £>t = [£ D + i VD }/^/2. This 
provides new ways to recombine the Majorana fermions. 
Of the various alternatives, consider in particular the 
Dirac fermions / = [£,c+i£,D\/ V% and g — [rjD — if]c]/V2, 
which we choose to yield real coefficients later on. Again 
suppressing the site index, this transformation may be 
summarised by 



C*-i[/ + / t +.g-.9 t ], 

£> = f[/ t -/-<?-s t ], 



(2.7) 



and its inverse 



f = ±[C + & + iD + iD^}, 
g=l[C-Cft + iD-iDl]. 



(2.8) 



As an aside, we note that the Jordan- Wigner transforma- 
tion applied to two copies of the Ising model does not by 
itself generate the correct commutation relations between 
pairs of spin operators a x taken one from each copy. To 
ensure these commutation relations one should in addi- 
tion introduce Klein factors. Since the Klein factors ul- 
timately have no effect on the equations we present, we 
omit them throughout this paper. 

For the doubled system, we are concerned with the 
transfer matrix products (suppressing the slice index) 
Hg Hp and VcVd- The value of the transformation 
Eq. (2/7) is that it reduces these products to the simple 
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forms 



H c Hd = exp(-2E^! <Mh + fi9i\) 
VcVd = cxp ^2 J^Zi 1 K> n ,i[g\f i+ i + ft + i9i] + B 



where the boundary term B is 



B 



(e ivNc + ^ ND )W M h + ft9M 



iirNc _ p inN D \t n i ft 



e^MJi+figM] 



(2.9) 



(2.10) 



This process of doubling the degrees of freedom and 
rewriting them locally as fermions, in order to remove 
terms which are not particle conserving, may be viewed 
as a local Bogoliubov transformation. 

The boundary term B contains the two boundary op- 
erators 



B 1 



e inN c _j_ g i7rJV D 



(2.11) 



These operators commute with the transfer matrix as a 
consequence of Z 2 symmetry: for a single system, say C, 
one can identify two invariant subspaccs, distinguished 
by the behaviour of vectors within the subspace under 
the operation Rc "which reverses the orientation of a com- 
plete row of spins B Specifically, 



Rc<J X jC R C = -VjC 



(2.12) 



for all j, and R c — 1. Introducing the corresponding 
operator Ro for the D system and assuming the total 
number of spins across the strip to be even, one finds that 
the boundary operators are simply — Rc±Rd- Since 
both Rc and Rjj commute with the transfer matrix, four 
invariant subspaces arise naturally from [Rc = ±1] (£> 
[Rd = ±1]. Using obvious notation, TcTjy may then be 
presented schematically in the block-diagonal form 



TcT, 



c-ld 



f ++ \ 
-+ 
+- 



V 







(2.13) 



/ 



Thus the Fock space associated with the C and D 
fermions can be divided into four subspaces according 
to the parity of Nc and Njj. In two of them, for which 
B~ = 0, the number of / and g fermions is conserved 
under the action of T. 



B. Network model interpretation 

The conservation of the Dirac fermions / and g under 
the action of the transfer matrix operator makes it pos- 
sible to go from a second-quantised description to a first- 
quantised form. Moreover, just as the second-quantised 



form has SO (2) symmetry.^ one finds that the first- 
quantised form may be interpreted as the transfer matrix 
for a scattering problem, because it fulfills the require- 
ments arising from unitarity of the scattering matrix. 
Specifically, the first-quantised form represents a network 
model, in which non-interacting / and g fermions propa- 
gate on directed links of a lattice. The fermions scatter at 
nodes, where two incoming links and two outgoing links 
meet. In this way, the nodes of the network model take 
the place of bonds in the Ising model. A correspondence 
of this type was set out first by Cho and FisherQ and 
subsequently refined by Cho,B who pointed out that the 
network model studied numerically in Ref. ^ is equiva- 
lent to an Ising model in which some exchange couplings 
are imaginary, while the RBIM itself is represented by 
a different network model. In this subsection we review 
these ideas. 

The identification of the first-quantised form of T 
makes use of a general equivalence between first- 
and second-quantised versions of linear transformations. 
Consider in a Hilbert space of dimension N a linear 
transformation of single-particle wavefunctions, repre- 
sented in a certain basis by an N x N matrix with el- 
ements (expG)y . Introducing in the same basis fermion 
creation and annihilation operators, a\ and o^, the 
second-quantised representation of this transformation is 
expfajGjjCHj]. To apply this equivalence to the transfer 
matrix T, let the N = 2M fermion annihilation opera- 
tors be [ai, • • • , a 2 Af] = [ft, • • • , /m, fll, ■ • • , 3m]- In the 
B~ = subspaces, the transfer matrix of the RBIM has 
the canonical form 



T C T D = ex.p[alGijaj], 



(2.14) 



and can be represented equivalently by the 2M x 2M 
matrix T, with elements Ty — (exp G)ij , as a transfor- 
mation of single-particle states. Thus the action of the 
operator T on a Slater determinant is replicated by the 
action of the matrix T on the orbitals entering the deter- 
minant. In the following we use notation for the matrix 
T corresponding to that introduced for the operator T: 
T n denotes the transfer matrix for the n-th slice of the 
system, T(k, I) indicates a product and T is shorthand 
for either. 

While knowledge of the single-particle form of T is 
enough by itself for efficient numerical calculations, phys- 
ical interpretation within this framework of the RBIM 
as a localisation problem depends on the fact that T is 
a pseudo-orthogonal matrix. In consequence, it can be 
viewed as the transfer matrix for a scattering problem 
in which flux is conserved. In order to see that this is 
indeed the case, consider the basic building blocks of 
the transfer matrix for one column of sites in the dou- 
bled Ising model. The two factors, HdHc an d VdVc, 
appearing in Eq. ( |2.2| ) each consist of products of M 
commuting operators. Every such operator represents 
a single bond of the Ising model and involves only one 
pair of / and g fermions. Schematically, a horizontal 



4 



bond gives rise to exp(— 2K*[g t / + pg]), which is re- 
placed in a first-quantised treatment by the 2x2 ma- 
trix h = expf— 2k*ct x ), while a vertical bond yields 
exp(2/c[<7 T / + yg\), which is replaced by v = exp(2K<j x ). 
To arrive at a scattering problem, the / fermions are re- 
garded (arbitrarily) as right-movers, and the g fermions 
as left-movers. Then the matrices h and v are transfer 
matrices for nodes of the network model. They relate 
flux amplitudes, Li„ and L out , to the amplitudes Ri n 
and R ut, appearing either side of a node as illustrated 
in Fig. |^. In algebraic terms, we have for horizontal bonds 



a) 




b) 




'-'out ±l in *~1n ±l out 

FIG. 2. Scattering nodes for: (a) horizontal and (b) vertical 
bonds. 



cosh 2k* — sinh2ft* 
-sinh2K* cosh 2k* 





and for vertical bonds the equation 

/ Rin \ ( cosh 2k sinh2«\ / L ou t 
\ Rout ) \ sinh2re cosh 2k I \ Li n 



(2.15) 



(2.16) 



Flux conservation follows from the relations o~ z b)o z = 
ft, -1 and cr z v^<J z = v~ l . 

The network model as a whole is illustrated in Fig. |[ It 
has the same structure as the U(l) network model, intro- 
duced to describe localisation in the context of the inte- 
ger quantum Hall effect .E^l Directed links form plaquettes, 
each with a definite sense of circulation, which is alter- 
nately clockwise and anticlockwise on successive squares. 
Disorder appears in the U(l) network model in the form 
of quenched random phases associated with links. By 
contrast, for the RBIM randomness enters only through 
the scattering parameters, 2k and 2k* , associated with 
nodes. An antiferromagnetic vertical bond leads to a 
negative node parameter, k. An antiferromagnetic hor- 
izontal bo nd, however, gives rise to a complex k*, since 
from Eq. (O) 



■ iw/2 , 



(2.17) 



generating an overall minus sign for h. The sign is ac- 
companied by a mi nus sign as a factor in the coefficient 
A 2 , defined in Eq. (fj). 

The form of this disorder determines the symmetry 
class to which this network model belongs hij-the classi- 
fication introduced by Altland and ZirnbauerEj. Specif- 
ically, Hamiltonians H belonging to class D have, in a 




FIG. 3. The network model. Flux propagates on links in 
the direction indicated by arrows. The transfer matrix relates 
flux amplitudes carried by links on the right to those on the 
left. Nodes arising from single rows of vertical bonds and 
horizontal bonds in the Ising model are indicated by V and 
H, respectively. Two particular nodes are labelled by h and 
v. Four sites of the Ising model are also shown with exchange 
interactions as dotted lines. 



suitable basis, the property that H* = —H, so that H 
is pure imaginary. Adapting this defining relation to a 
network model, one supposes that propagation on the 
network is generated by a time-evolution operator for 
unit time-step, exp(iH). For class D, this evolution op- 
erator is real, so that scattering phase factors may take 
only the values ±1, as is indeed the case for the RBIM. 
In detail, a single antiferromagnetic bond (either hori- 
zontal or vertical) introduces phases of it for propagation 
around both the anticlockwise plaquettes that meet at 
the corresponding node, compared to the phases in the 
purely ferromagnetic model. Other choices of random- 
ness belonging to the same symmetry class are of course 
possible. Cho and Fisheru investigated a model in which 
the transfer matrices at all nodes are of the type given 
in Eq. (2.16), n wj.th randomness in the sign of K, while 
other authorsErEJ have studied a model in which scatter- 
ing phase factors of ±1 are associated independently with 
links rather than nodes. Strikingly, each of these differ- 
ent choices leads to vexiz-different localisation properties 
in the network model.Qo 

Combining the 2x2 transfer matrices, h or v, for each 
node, one arrives at the 2M x 2M transfer matrix T for 
the system as a whole. Flux conservation guarantees that 
T may be factorised as 



T = 




cosh(eL) sinh(eL) 
sinh(eL) cosh(eL) 



v R J 

(2.18) 



where components in the basis are ordered so that the 
amplitudes for propagation in one direction constitute 
the first M entries of the vectors on which T acts, and 
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those for propagation in the opposite direction make up 
the remaining M entries. Here, the M x M matrices, 
Wl, Wr, Vl and Vr are for a general localisation prob- 
lem unitary matrices, and for the Ising model orthogonal 
matrices, since in that case every element of the trans- 
fer matrix is real. The M x M matrix e is r eal, p ositive 
and diagonal. It is convenient to rewrite Eq. ( [2.18 ) in the 
form 



l ( W L -W L 








-cL 



As a starting point, consider the polar decomposition 
of the transfer matrix for the doubled Ising model, which 
takes the form 

2 M 

T c f D = J2 \ L *c) ® \L jD )e^ +x ^ L {R jD \ <g> {Ric\ . 

(2.22) 

y R \ Here, {\L iC ) <E> \L jD )} and {(Rjd\ <8> (Ric\} are two com- 

I ] plete, orthonormal sets of many-particles states for the C 
-Wr Vr J and D fermions, which in general are not bi-orthogonal. 

(2.19) The factors e XiL are the singular values of the transfer 



where the diagonal elements of exp(±eL) are the singu- 
lar values of T. For a random system of length L, the 
exponents eL are O(L), with sample-to-sample fluctua- 
tions which are ©(L 1 / 2 ). From Oseledec's theorem, the 
average e tends to a limit, diag(ei, e%, ■ ■ ■ , £m), for large 
L, where ei < £2 < • • ■ £m are the Lyapunov exponents 
characterising the network model. 

It is useful also to express Eq. (2.19) in second- 
quantised notation. Writing the left and right orthogonal 
matrices in terms of the Hermitian 2M x 2M matrices 
Al and Ar, defined by 



exp[-iA L ] = -±5 



exp[iA fl ] = 



W L -W L 

wl V R \ 
-Wr Vr) 



(2.20) 



the transfer matrix for the doubled Ising model takes the 
form (within the subspaces with B~ = 0) 



f c f D = exp[-ia\ AujOtj] x exp[a\ai[a z ® eL] u ] 
x exp[ia\ARijaj] . 



(2.21) 



C. Lyapunov exponent spectrum 



In this subsection we discuss some aspects of the map- 
ping between the RBIM and the network model that 
have not been considered in previous work. These stem 
from the fact that, under the Jordan- Wigner transforma- 
tion, different boundary conditions arise in T according 
to t he pa rity of the fermion numbers Nc and No (see 
Eq. ( 2.1 0| ) ) . Full information on sectors of both parities 
is contained in the results of network model ca lcula tions 

for the subspaces denoted ++ and in Eq. ( 2.13 ). To 

make use of this information it is necessary establish how 
the Lyapunov exponents of the spin transfer matrix are 
related to those of the network model. A crucial step 
is to be able to identify the parity of left and right vec- 
tors of T when these are written in terms of the / and g 
fermions. We show here how this may be done. 



matrix for a single copy of the spin system, and the lim- 
iting values of for large L are the Lyapunov exponents 
characterising the spin system. For economy, we use the 
same symbol to denote both the disorder-dependent X{ 
at finite L and its limiting value as L — > oo. Since we are 
concerned with the largest few singular values, we adopt 
the ordering Ai > A2 > • • • > A 2 j 



Comparing Eq. (2.21) with Eq. (2.22), one sees that the 
values taken by exp(aj a.i[a z ® eL]a) for a\ai = or 1 
generate the possible values of eS Xi+Xj ^ L . In particular, 
ignoring for the moment questions connected with parity, 
the largest of the Lyapunov exponents for the doubled 
Ising model is obtained by setting a\ai — 1 for 1 < i < 
M and a\ai = for M + 1 < i < 2M. The associated 
right vector is 



I*) = \Ru 



\Ru 



M 

=n^[^/i+^}]i vac ) 



(2.23) 



where |vac) is the vacuum for / and g fermions, and for 
simplicity we have omitted the subscript R on W and V. 
The state |\&) satisfies for all i the equations 



Wg/t + V i3 g\ 







and 



Mfj - v ij9j ] 1*) = 0. 



(2.24) 



(2.25) 



Let P = l[W T -V] and Q = \[W T + V]. Taking 



the difference between Eq. ( p. 24 ) and Eq. ( p. 25 ) yields 



7i = for all i, where the fermion creation opera- 
tors 7J are defined by 



= P ij C j + > 



lijC r 



(2.26) 



(Of course, similar expressions for t he D syste m m ay be 
obtained from the sum of Eq. ( |2~24| ) and Eq. In 
this way we find that the right vector associated with the 
largest possible singular value of the spin transfer matrix 



M 



\Ric 



INio>, 



(2.27) 



G 



where |0) is the vacuum for the 7-fermions. More gen- 
erally, we can obtain all the right vectors as follows. 
First, in the factor exp(a\ai[(j z ® eL]a) from Eq. ( 2.21 ), 
for each i in the range 1 < i < M either: (a) set 
a\on — 1 and a\, M oti+M — 0; or (b) set a\a.i — and 
a|, M a>i+M = 1- The corresponding right vector \R) sat- 
isfies for (a) i}\R) = and for (b) ji\R) = 0. The as- 
sociated Lyapunov exponents for the (undoubled) Ising 
model are 



M 



(2.28) 



where jjji = 1 or for (a) and (b) respectively. 

As a further step in the discussion, it is necessary to 
distinguish between the two sectors with even and odd 
parity for the fermion numbers Nc and Nr>. Except in 
strip geometry (k^m = in Eq. (2.10)), different bound- 



ary conditions are imposed on the network model for each 
sector, and so each sector has its own set of Lyapunov 
exponents, e, and matrices, W and V. We indicate quan- 
tities calculated using boundary conditions appropriate 
for even and odd parity sectors with plus and minus signs 
respectively: e^, and V^. Introducing the number 
operator for 7 fermions, N 7 — Y^iLi^Hlh ^ ^ s straight- 
forward to see that, in general, either exp(j7r./V c ) = 
exp(z7riV 7 ) or exp(i7riV c ) = — exp(i7riV 7 ), but to deter- 
mine which of these holds in a particular instance re- 
quires explicit (numerical) calculation. To this end, we 
consider (restricting ourselves for simplicity to even M) 
the scalar product of |^) (see Eq. (2.23)), for which we 
know that e mNri \ty) = +| V E'), with a reference state, |ref), 
chosen in order that exp(i7rA r c )|rcf) = +|ref). The result 
(ref|W) 7^ will indicate e l7rNc — e ITrJV T, while (barring 
accidental orthogonality) the result (ref | = implies 
exp(z7rA r c ) = — exp(z7riV 7 ). A suitable choice for |ref) is 
the state 



M 



|rcf) =[](//+ sl)|vac), 



(2.29) 



which satisfies 7Vc|ref) = -/Vo|ref) = M|ref) and hence 
also exp(z7rA r c )|ref) = +|ref). The scalar product is 



(ref|tf> = 2- M / 2 det(iy T + V) 

= 2- M / 2 det(iy)det(l + WV) . 



(2.30) 



The only factor on the right side of this expression 
which may be zero is dct(l + WV). It turns out that 
X = det(WV), which takes the values \ = ±1, is a con- 
venient indicator: barring accidental degeneracies in the 
spectrum of W V , det ( 1 + W V) = if and only if x = - 1 • 
The proof of this statement is as follows. One has 



det(l + WV) = ]J(1+ Pi ) 



(2.31) 



where pi are the eigenvalues of the O(M) matrix WV. 
These occur as complex conjugate pairs, pi and p*, and 
possibly also as real pairs, 1 and —1, of which there will 
be at most one in the absence of degeneracy. If there is 
one such real pair, \ = — 1 and det(l + JW) = 0; if there 
is none, x = +1 and det(l + WV) ^ 0. 

We now apply these results to obtain expressions for 
the Lyapunov exponents of the Ising model transfer ma- 
trix in terms of those of the network model. For sim- 
plicity of presentation we make use of a property which 
appears to hold generally and is certainly true for the 
model studied in Sec|V|, the ±J-RBIM with p < 0.5. 
In this system, x + = +1 always, and hal f of t he Lya- 
punov exponents Ai are obtained from Eq. (2.28) by set- 



ting e = e + and taking A^ 7 even. The remaining expo- 
nents result from setting e = e - , accompanied by even 
A 7 if x~ — +lj an d by odd iV 7 if x~ = — 1. Since we 
are concerned in the following only with x~ > we write it 
below simply as x- 



Using the expression for the exponents, Eq. (2.28), we 
find the following rules for the case x = 1 



(2.32) 



Ai = 


1 + 

2 Zjj=l fc i 




A 2 = 


2 Zjj=i H t i 




A.3 = 


1 sr^M - 

2 2_,j=l £ i € 2 




A 4 = 


1 T M + + _ 
2 2^i=l H e l 





For the case x = — 1> we have instead 



Ai = 


1 s-^M + 

2 2^i=l fc i 


A 2 = 


ly M e - 


A3 = 


2 2-<i=l H 


A 4 = 


1 ^« ,+ 

2 Z^i=l fc i 



(2.33) 



where the order of A3 and A4 has to be decided numeri- 
cally. 

It is interesting to note a consequence that follows from 
the importance of x> an d which is probably characteris- 
tic of localisation problems in class D. It arises if x can 
change sign as a continuous parameter, such as temper- 
ature in the Ising model, is varied. Since the two sub- 
spaces of WV £ O(M) in which x = +1 an d X = — 1; 
respectively, are disconnected, a change in the sign of x 
is accompanied by the vanishing of . This process is 
a form of level crossing, as illustrated in Fig.^. In the 
RBIM it occurs for large M at the Curie point, as dis- 



cussed in Sec. IV 



This distinction between phases with either sign for x 
is the analogue for the RBIM in cylindrical geometry of a 
topological classification introduced for two-dimensional 
systems from class D in Rcf. p9| a nd for one-dimensional, 
single channel systems in Ref. [42|. In particular, such one- 
dimensional systems may have two phases: in one phase 
a long sample supports a zero-energy state at each of its 
ends, and in the other it does not. Turning to the network 
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FIG. 4. A sign change of x as a function of a parameter x 
is accompanied by the smallest Lyapunov exponent reaching 
zero. This may be regarded as a form of level crossing, as 
illustrated. 



model for large L, we note that the combinations V^W^ 
and WrVr are the reflection matrices from either end of 
the system. A closed sample may be constructed in an 
obvious way, by joining outgoing links to ingoing links in 
pairs at each end of the system. For a network model, 
a stationary state has the status of a zero energy state, 
and stationary states will exist at the ends of the closed 
sample if the reflection matrices for the corresponding 
open system h ave 1 as an eigenvalue. From the discussion 
following Eq. one sees that this is the case if % = 

— 1 but not if X = +1- 



III. CALCULATIONAL METHODS 

A. Numerical procedure 

Numerical methods suitable for studying random 
transfer matrix products in general are well-established 
and described, for example, in Refs. |4^-^. It has been 
recognised recently,Ea however, that these methods may 
develop an instability to rounding errors and must be 
modified when applied to systems in symmetry class D. 
Specifically, the modifications are required if the smallest 
positive Lyapunov exponent approaches zero on a scale 
set by the spacing between other exponents, which hap- 
pens in the RBIM at the Curie point, as described in 
the Sec. II C and Sec. IV. We summarise the established 



algorithm and review the modification required in this 
subsection. 

First, we define some notation. Consider a network 
model of width 2M links and len gth £ , with a transfer 



matrix of the form given in Eq. (OS). Let x (n), for 



k = 1,2,- •• , 2M and n fixed, be orthonormal column 
vectors, each of 2M components, written in the same 
basis as this transfer matrix. These vectors arc generated 
by a sequence of operations designed to ensure that x fc (L) 




converges for large L to the k-th column of the matrix 



(3.1) 



appearing in the polar decomposition, Eq. ( p,18j ). 

The conventional choiceoc3 for these operations is as 
follows. Pick x fe (0) arbitrarily. With n = 0, let 

y fe =T{L-n-s 1 L-n)yL k (n) 1 (3.2) 

and perform Gram-Schmidt orthonormalisation, follow- 
ing 

k-l 

z k = y k _ £ ([x i(„ + s) f . y fc )x> + s) (3 .3) 



i=l 



and 



x fe (n - 



(3.4) 



The process is repeated with n = s,2s . . . L — s. The 
Lyapunov exponents are then the mean growth rates 



,1 



Eft 



ln|z M+1 



-fci 



,1 



ln|z M+fc | 



(3.5) 



for k — 1 . . . M, where the average is over successive or- 
thonormalisation steps. The interval s is taken for com- 
putational efficiency to be as large as is possible without 
rounding errors significantly affecting the orthogonalisa- 
tion. 

The rate of approach with i ncr easing L of the vectors 
x'(L) to the columns of Eq. (3.1) is determined by the 
spacing between successive Lyapunov exponents. So also 
are the deviati ons at large L of these vectors from the 
columns of Eq. ( |3.l|) . Such deviations are induced by nu- 
merical noise and generate errors in the calculated values 
of Lyapunov exponents. For systems in symmetry class 
D, the value of the smallest positive Lyapunov exponent, 
ei, may approach zero. If it does, the vectors x M (L) and 
x M+1 (L) are unusually susceptible to r ound ing errors, as 
is the value of t\ determined from Eq. (|3.5[ ). We demon- 
strate in Appendix |A| that the error decreases with de- 
creasing noise amplitude, ct, only as |ln(cr)| _1 . Because 
of this, a modification must be found that stabilises the 
algorithm. 

Following Ref. |4l], we adapt the Gram-Schmidt or- 
thonormalisation to enforce the 2x2 block structure evi- 
dent in Eq. (3.1). Denoting the j-th component of x fc (n ) 
by (n), and similarly for y fc and z fe , we replace Eq. ( |3.3| ) 
for 1 < k < M by 



fc-l M 



i=l ' 1=1 

if 1 < j < M, and by 

fc-l 2M 



(5>K" + s)yf)x i j (n + s) (3.6) 



*i=Wi-E( E x\{n + s)yi?)x){n + s) (3.7) 



i=l l=M+l 
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if M + 1 < j < 2M. Similarly, we replace Eq. ( J3T4| ) by 



fe|2 



1/2 



(3.8) 



if 1 < j < M, and by 



2M 



x k j (n+s) = z k j /[ Yl 



_fc|2 



1/2 



(3.9) 



if M + 1 < j < 2M. Lyapunov exponents are deter- 
mined as before from Eq. ( |3.5D , and now remain stable to 
rounding errors even if t\ — ► 0. 



B. Self-averaging quantities 

We wish to calculate for the Ising model the free en- 
ergy, spin correlation functions, and correlations of dis- 
order operators. In the presence of bond randomness 
these all exhibit sample-to-sample fluctuations, but the 
free energy density and typical decay lengths appearing 
in correlations functions are self-averaging. In this sub- 
section we describe how such self-averaging quantities can 
be obtained from the Lyapunov exponent spectrum of the 
network model. The calculation of correlation functions 



themselves is discussed in Sec. [II C 



We start from the polar decomposition of the transfer 
matrix for an (undoubled) Ising m odel of width M and 
length L, which (in analogy to Eq. ( 2.22| )) is 



T 



E 



|^)e A;i (i?i 



Defining the reduced free energy per site as 



F 



lim ln(Z)/LM 



(3.10) 



(3.11) 



and using Eq. (gjj), Eq. (|2.32j) and Eq. fl2.33|) , we have by 
standard arguments 



F = - lim \ 



In \A\ 



M 

— Y 

2M ^ 



(3.12) 



Turning our attention to typical decay lengths, we note 
first that, viewing the network model as a localisation 
problem, its smallest positive Lyapunov exponent defines 
a localisation length £ through £ = e7 • In a localised 
phase £ has a finite limit, the bulk localisation length, 
as M — > oo, while at a mobility edge one expects that £ 
diverges with M and that a universal scaling amplitude, 
a, is defined by the limiting value of Mt\ for M — > oo. 
An unusual feature of localisation problems in symmetry 
class D is that one may have a — 0; from the discussion 



of Sec. 1IC and results presented in Sec. IV, this occurs 



in the RBIM in the sector with odd parity. 



For the Ising model, the typical correlation length £ aa 
appearing in the spin-spin correlation function may be 
extracted as follows. This correlator, for two spins with 
(in the notation of Fig. |l|) the same coordinates in the 
vertical direction and separation n — m in the horizontal 
direction, is 



Tr[<rfT(l,n)(jff(n+l,L)] 
Tr[T(l,L)] 



(3.13) 



Recalling that af has non-zero matrix elements only be- 
tween states with opposite parity, and taking L — * oo, 
^uu is defined and expressed in terms of the Lyapunov 
exponents for the spin transfer matrix by 



Ca=~ lim -ln(«(0)<H)) = A 1 



A 2 



(3.14) 



When writing £ CT(T in terms of the network model Lya- 
punov exponents, it is useful to introduce a lengthscale 
£id which characterises the sensitivity of the network 
model to changes in boundary conditions, and is defined 
by 



£id " 2 5Z[ e *~ 



(3.15) 



We expect insensitivity to boundary conditions except at 
the critical point, and anticipate that ~ exp(— M/£) 
for large M . In regions of the RBIM phase diagram for 
which x = 1 (corresponding, as we argue, to the param- 
agnet), we have from Eq. (2.32) 



1 

aa 



•S3 



(3.16) 



so that asymptotically the localisation length, £, and spin 
correlation length, £ CTCr , are equal. By contrast, in regions 
of the phase diagram for which % = — 1 (corresponding to 
the ferromagnet) we have £ CTCr = £id- This large length- 
scale here characterises the decay of spin correlations in 
a quasi-one dimensional sample within the ordered phase 
of the two-dimensional system. Such decay is governed 
by rare domain-wall excitations that cross the width of 
the sample. Because £ CT(T is large when \ = — L h is 
useful also to exa mine t he inverse lengthscale governing 
corrections to Eq. (3.14), which is Ai — A3. For x = — 1 



Ai -A 3 



ID ' 



(3.17) 



so that, as M 



e 2 gives the typical decay rate 



of the connected part of the spin correlation function in 
the ordered phase. 

In a similar way, one can obtain the typical cor- 
relation loBgth for the disorder operators fi r of Kadanoff 
and CevaJ£3 These operators are defined at points r which 
lie at the centres of plaquettes in the Ising model. The 
two-point correlation function (^ofi r ) is defined by con- 
sidering a modified system in which exchange interactions 
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crossed by a path on the dual lattice between and r have 
their sign changed. Then 



(HOfJbr) = 



~Z 



(3.18) 



is the ratio of the partition function Z' for the modified 
system to that of the original system, and 



i 



lim - ln({/ioMr)) 

r — >oo t 



(3.19) 



Because the different boundary conditions imposed on 
the network model in sectors of even and odd parity con- 
stitute an (infinite) line of such modified bonds, may 
be expressed in terms of e + and e~. Moreover, in the fer- 
romagnetic phase (x = —1), £ w is the reduced interfacial 
tension between domains of opposite magnetisation. To 
make this explicit, let F p and F a be reduce d free energies 
per site, calculated from the definition Eq. ( 3.12|) for sys- 
tems in cylindrical geometry with, respectively, periodic 
(°m+i — a i) an( i antiperiodic (cr^ +1 = —erf) boundary 
conditions on spins imposed around the cylinder. Then 

C = M(F P - F a ) . (3.20) 

In this phase, we find using the ideas of SecETO that 



£uu — e i 



(3.21) 



while in the paramagnetic phase (\ = +1) we obtain 
= £id, so the decay length diverges with M. As 
one might expect, the behaviour of £ w in each phase is 
similar to that of ^ aa in the dual phase. 



C. Correlation functions 

Calculation of the full form of correlation functions 
is more involved than that of the typical decay lengths 
since, of course, the results cannot be expressed solely in 
terms of Lyapunov exponents. Nevertheless, it turns out 
that even powers of correlation functions may be deter- 
mined straightforwardly.!] In the most important exam- 
ple of the second power, one requires the product of two 
equivalent correlation functions, evaluated for each of the 
two copies of the Ising model that are combined in the 
network model. In the case of the square of the two- point 
correlation function of disorder operators, Eq. (3.18), this 
means that the same modification of bonds is introduced 
in both copies of the Ising model, so that (Z') 2 is de- 
termined from a network model with a specific set of 
modified nodes. In the case of the square of the spin- 
spin correlation function, one can take a similar route by 
expressing this in terms of a disorder correlator in a dual 
system. Alternatively, one can write the product of two 
copies of a spin operator in terms of / and g fermions, 
as we describe below. By either route, one arrives ulti- 
mately at the same result: the square of the spin-spin 



correlation function is given by the ratio of the square of 
a partition function calculated from a modified network 
model to same quantity calculated from an unmodified 
model. By contrast, odd powers of correlation functions, 
including the first power, appear to be much harder to 
evaluate, leaving the sign of the correlation function un- 
determined: we summarise the difficulties that arise at 
the end of this subsection. 

To obtain t he s quared spin-spin correlation function 
following Eq. (3.13), we must evaluate products involving 
the transfer matrix for the doubled Isin g m odel and also 
factors of the form a x c a x D . From Eqs. ( |2.4| ) and (2.7) we 
have 

= exp (in YtfUi + alfi] + i*f]f 3 

1=1 

(3.22) 



In the spirit of Sec. II B we translate this into first- 
quantised form. Each operator exp(in[fjgi + gjfi]) is 
represented by a 2 x 2 matrix, exp(iira x ) = — 1. As a 
result, on one slice of the network model phase factors of 
— 1 are associated with each of the right and left going 
links having coordinate i in the range 1 < i < j — 1. 
In addition, the operator expiiirf -fj) is represented by 
a similar phase factor associated with the j-th right- 
going link. These phase factors are illustrated schemat- 
ically in Fig. (p^,), using as an example the combina- 
tion <J x c a x D a% c a %Di wn i cn arises in the calculation of 
(af(n)a x (m)) 2 on setting n = m, i = 1 and j = 4. 
Such link phases can equally be attributed to nodes rep- 
resenting vertical bonds of the Ising model, as indicated 
in Fig. dab). Viewed in this way, the insertion of spin 
operators into the transfer matrix product is represented 
by a change in node parameters n n j — > n n j + in/2 for 
1 < i < j ' — 1. In turn, this is equivalent to a change 
in sign for the corresponding dual bond strengths, as it 
should be since the spin correlation function can be eval- 
uated as a disorder correlator in the dual model. 

Implementing this approach in numerical calculations, 
we determine the singular values of the transfer matrix 
for modified and unmodified network models of length 
L, of course using the same realisation of disorder for 
both. From these we calculate the largest singular value 
of the transfer matrix for the doubled spin system, which 
we denote by exp(2A' 1 L) in the modified case, and by 
exp(2Ai_L) in the unmodified case, following the notation 
of Sec. [II Q. For large L 



(annK(m)) 2 =exp(2[A' 1 -A 1 ]L). 



(3.23) 



In practice, the combination [A^ — Ai]Z approaches a finite 
limiting value rather quickly with increasing L. Conve- 
niently, it is not necessary to evaluate the scalar products 
of the form |(Li|i?i)| 2 whi ch appear in the numerator and 
denominator of Eq. ( 3.13| ), because for large L these are 
the same in the modified and unmodified systems, and 
therefore cancel. 
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probability distribution 

P(Jy) = (1 - p)S(Jij - J)+p6{J l3 + J) . 



(4.1) 



(a) (b) 

FIG. 5. Schematic representation of the effect of modifying 
the system by inserting spin operators in one slice: (a) link 
phases associated with a pair of spin operators; (b) equivalent 
node phases. 



As mentioned above, calculation of the unsquared spin- 
spin correlation function presents greater practical prob- 
lems. A route is cl ear in principle: one can use the dis- 
cussion of Sec. |ll C| to construct the transfer matrix for 
the undoubled Ising model, via its polar decomposition, 
in terms of Slater determinants of the 7 fermions; and 
one can also express spin operators in this Ising model 
in terms of the creation and annihilation operators for 
these fermions. Difficulties then arise from the fact that 



the matrices Wl and Wr appearing in Eq. (2.19) are un- 
related in the presence of disorder, as also are Vl and 
Vr. In consequence, one has to deal with two sets of 7 
fermions: jl and 77?. Put briefly, we find (as in Eq. (4.7) 
of Schultz, Mattis and LiebO) that (of (n)<Jj (n)) can be 
written as an expectation value of a product of 2\i — j\ 
fcrmion operators, which can be evaluated using Wick's 
theorem. However, in the disordered system it is not 
possible to reduce this expectation value to a single de- 
terminant (as in Eq. (4.13) of Ref||). Without such a re- 
duction, the computational effort required to determine 
(of (n)oj(n)) seems prohibitive for large \i — j\. 



IV. NUMERICAL RESULTS FOR THE ±J-RBIM 



A. Introduction 



with < p < 1 and J positive; we set J — 1 in the 
following. 

The phase diagram of the model, as a function of 
temperature T and the concentration p of antiferro- 
magnetic bonds, is shown in Fig.^|, with renDrjmalisa- 
tion group (RG) scaling flow superimposed.E3~o The 
pure system {p — 0) has a transition between ferro- 
magnetic and paramagnetic phases at a Curie temper- 
ature Tq — 2[ln(l + \/2)] _1 - As antiferromagnetic bonds 
are introduced the Curie temperature is depressed, and 
the ferromagnetic phase is destroyed altogether above a 
threshold concentration p c . A, cairve in the p — T plane 
known as the Nishimori lineEfltll (NL) plays an impor- 
tant role in the discussion of scaling flow. It is defined 
for the ± J RBIM by the equation exp(2/3J) = (l-p)/p. 
On this line the RBIM has an additional gauge symme- 
try, because of which the internal energy is analytic and 
ensemble-averaged spin-spin correlations obey the equal- 
ities [(of (n)oJ(m)) 2fe_1 ] = [(of (n)oj (m)) 2k ] for integer 
k. The NL cuts the phase boundary separating the ferro- 
magnet from the paramagnet at a point C, the Nishimori 
point, with coordinates p c , Tjy. This point is particularly 
interesting as an example of a disorder-dominated mul- 
ticritial point. One of the two scaling flow axes in its 
vicinity lies along the-.NL, while the other coincides with 
the phase boundaryjlj as indicated in Fig. ^. Scaling flow 
on the critical manifold for p < p c runs from the Nishi- 
mori point towards the critical fixed point of the pure 
system, at which disorder is marginally irrelevant. The 
phase boundary on the okhc£.side of the Nishimori point 
is believed to be verticaOlLZl in the T — p plane, and on 
it the scaling flow runs from the Nishimori point towards 
a zero-temperature critical point. Finally, the phase di- 
agram for p > 1/2 can be obtained from that shown for 
p < 1/2 by reflection in the line p = 1/2, using a gauge 
transformation which maps p to 1—p and the ferromag- 
netically ordered phase to an antiferromagnct. 

Despite the considerable effort which has been invested 
in studies of the RBIM, some aspsects of its behaviour are 
not yet well-characterised. In the following, we present a 
high-accuracy determination of the position of the phase 
boundary and of critical properties at the Nishimori 
point. 



In this section we present results obtained using the 
mapping from the Ising model to the network model as 
a way of studying the ± J RBIM. Previous work of this 
type has been described by Chocl, but without the advan- 
tages of the numerical algorithm or the detailed relation 
between the network model and statistical mechanical 
quantities that we have discussed in Sec. III. The ±J 



B. Method 



RBIM, defined on a square lattice, has nearest-neighbour 
exchange couplings drawn independently from the 



We use the numerical method described in Sec. Ill A 



to 

calculate the Lyapunov exponents of the network model 
associated with the RBIM, studying two copies of the 
system for each disorder realisation, with boundary con- 
ditions appropriate for fcrmion numbers of each parity. 
In the spirit of Sec. [II B we use the smallest positive 
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P c p 0.5 

FIG. 6. Phase diagram of the ±J-RBIM with superim- 
posed RG scaling flow. 



For this purpose the quantity %, introduced in Sec.IIC, 
is very useful and we substantiate our claim that (in the 
thermodynamic limit) the sign of x indicates which phase 
the Ising model is in. 



2.0 



1.0 



0.0 



T=0.68+4.67p 

.---16' 



0.00 



0.10 

p 



0.20 



FIG. 7. The location of the phase boundary determined 
from numerical calculations. Data obtained on the line 
T = 0.68 + 4.67p are presented in Fig.| 



exponent calculated for the network model with peri- 
odic boundary conditions to define a characteristic in- 
verse lengthscale, and analyse the finite-size scaling be- 
haviour of M/£ = Mci as a function of system width 
M. In addition, we determine the interfacial tension, 
and study its size-dependence. We also calculate 
the disorder-averaged square of the spin-spin correlation 
function, [(af(n)a^(n)} 2 }, for spins lying in the sam e slice 



of the system, using the approach described in Sec. [Ill C . 

For most of the results presented, we study system 
widths in the range from M = 8 to M = 256 spins, 
and system lengths of L = 5 x 10 5 spins. Realisation- 
dependent fluctuations in self-averaging quantities de- 
crease as L" 1 / 2 and in some cases increase with M. As 
an example, using L = 5 x 10 5 the value of ej~ at the 
Nishimori point is obtained with an accuracy of 1% for 
M = 16 and 2% for M = 64. Some calculations require 
higher precision. In particular, the high-resolution stud- 
ies of the interfacial tension close to the Nishimori point, 
presented in Sec. IV 5 , and of scaling on the phase bound- 
ary, presented in Sec. IV E, use systems of length up to 
L = 2 x 10 s , restricting accessible widths to M < 32. 



C. Location of phase boundary 

In this subsection we describe the determination of 
the form of the boundary between the ferromagnetic 
and paramagnetic phases. We also discuss the nature of 
finite-size effects in different parts of the phase diagram. 



Our results for the position of the phase boundary are 
shown in FigJ^ and in Table |. Points on this phase 
boundary are found from a finite-size scaling analysis of 
the variation of M/£ along lines that intersect it; the 
slopes of these lines in the p — T plane are chosen to 
avoid crossing the boundary at small angles. Represen- 
tative data, calculated on the line T = 0.68 + 4.67p, are 
shown in Fig.^; they have two features that can be used 
to identify the boundary. First, the curves of M/£ for 
two successive values of M have an intersection point, 
and with increasing M these intersection points approach 
the boundary from the small-p side. Second, for each M , 
there is a value of p at which £ diverges, or equivalently 
(-i = 0. With increasing M, these points approach the 
boundary from the large-p side. We obtain consistent 
results using the two methods. 

A test of these calculations follows from the fact that 
the tangent to the ferromagnetic-paramagneticJpound- 
ary at the pure critical point is known exactlyca to be 
dTp/dp\ p= o = —7.2821.... From a linear approximation 
at p = 0.005 we find dT p /dp\ p=0 = -7.32 ± 0.06, in good 
agreement with this. Our values for T p are also compat- 
ible with those given in Ref. |2(| 

It is evident from the data shown in Fig.^J, and its 
equivalent for other values of p and T, that 6^=0 along 
a line in the phase diagram which approaches the phase 
boundary for large M, but is displaced from it into the 
paramagneti c ph ase for finite M. From the discussion 
given in Sec.IIC, we expect x t° change sign on this 
same line, being for large M positive in the paramag- 
netic phase and negative in the ferromagnetic phase. The 
data shown in Fig.^| demonstrates that this is so; Fig.^ 
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0.102 



FIG. 8. Values of 2M/£ calculated crossing the phase 
boundary along the line T = 0.68 + 4.67p. 



p 


T P 


P 


T P 


0.005 


2.2325 ± 0.0003 


0.0903 ± 0.0002 


1.458 


0.02 


2.120 ± 0.001 


0.0951 ± 0.0005 


1.379 


0.05 


1.875 ± 0.001 


0.1000 ± 0.0005 


1.294 


0.06 


1.783 ± 0.002 


0.1035 ± 0.0011 


1.224 


0.07 


1.688 ± 0.002 


0.1055 ± 0.0011 


1.173 


0.08 


1.580 ± 0.002 


0.1080 ± 0.0021 


1.095 


0.0852 


1.523 ± 0.002 


0.1090 ± 0.0021 


1.019 



TABLE I. Location of the phase boundary. 



also shows that the finite-size shift in the position of the 
phase boundary is very large in the portion of the phase 
diagram lying below the NL. It seems possible that these 
finite-size effects may provide an alternative explanation 
of data which have been-ipterpretecO'O as evidence for a 
random antiphase states lying in this region of the phase 
diagram; and it seems likely that they are responsible for 
non-monotonic temperature-dependence of Lyapunov ex- 
ponents for the RBIM, reported at p > p c in Ref. Eq. 



D. Nishimori Line 

In this subsection we examine critical behaviour near 
the multicritical point Cj-asJt is approached along the 
Nishimori line. The factst£fE3 that C is known to lie on 
the NL, and that the NL coincides with one of the scal- 
ing flow axes at C, both greatly help the analysis. We 
obtain consistent, high-accuracy estimates of the coordi- 
nate p c and the exponent v using three separate analyses 
of the finite-size scaling of £ and also from a study of the 
interfacial tension. 



0.4 
0.00 



0.05 



0.10 
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0.15 



FIG. 9. The sign of x f° r a system of width M = 16 as a 
function of position in the T — p plane. Open squares indi- 
cate x = +1 an d filled squares \ ■ = — 1. The NL and phase 
boundary are also shown. 



An overview of the variation of M/£ along the NL is 
given in Fig.|l0|. We apply finite-size scaling ideas to the 
data in the following different ways. Two of them are 
similar 



to the methods used in Sec. IV C to locate the 
phase boundary: first, curves of M/£ for two successive 
values of M cross, and we focus on these crossing points 
for increasing M; second, for each M there is a point on 
the NL at which M/£ = 0, and we study the position of 
these points as a function of M. Third, we can collapse 
data for different M and from the whole critical region 
onto a single curve. 

Turning to the first of these, we concentrate on the 
top left of Fig.jlO], where data sets intersect roughly at 
one point. Behaviour in this region is shown on a larger 
scale in Fig|ll]. From an extrapolation of the intersection 
points to large M we find p c = 0.1093 ± 0.0002. We 
also obtain a limiting value at the intersection point of 
M/i = 1.58 ± 0.01 as M -> oo. The value of v may 
be found from the scaling with M of the gradients of 
curves at the intersection points; a similar analysis can 
also be made for the interfacial tension and we present 
both together, towards the end of this subsection. 

Taking a second approach to the data, the points pm 
on the NL at which M/£ = are determined for 8 < 
M < 256 as shown in Figjl^, where we take advantage of 
the fact that, for fixed M, the combination xMj^ varies 
smoothly through zero as a function of position along the 
NL. One expects the finite-size shift pm —p c to vary with 
M as (pm — Pc) oc M~ x l v and we show the dependence 
of pm — p c on M in Fig. [13], using a double logarithmic 
scale for various choices of p c . With the correct choice for 
Pc this data should fall onto a straight line of slope —v. 
By this method we find p c = 0.1093 and v = 1.49 ± 0.05. 

A third treatment of the data for M/£ is provided by 
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FIG. 10. Variation of 2M/£ along the NL, with position 
parameterised by p. 
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FIG. 12. Variation of the combination 2\M/£ along the 



NL. 




0.1080 0.1090 0.1100 0.1110 



FIG. 11. Variation of 2M/£ along the NL, close to the 
Nishimori point. 
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FIG. 13. ln(j>M — Pc) as a function of ln(M), for differ- 
ent estimates of p c . The straightest line is obtained with 
p c = 0.1093 and has slope -1.49. 



attempting to collapse all points from the critical region 
of Fig. [To] onto a single curve, plotting M/£ as a function 
of (p— p c )M 1 ' v . In principle, both p c and v may be taken 
as fitting parameters, but we find that p c is more accu- 
rately determined using the methods described earlier. 
We therefore set p c — 0.1093 and vary only the value 
of v. We find the best collapse, shown in FigJTJ, tak- 
ing v — 1.50. Visibly worse collapse results from using 
v = 1.40, as shown in Fig.|l5]; by such comparisons we 
find v — 1.50 ± 0.10, confirming the result derived from 
Fig.[ll| but not improving on it in accuracy. 

Finally as a way to check the conclusions we have 
reached from finite-size scaling of M/£, and in order to 
make a direct coniBarison with recent work by Honecker, 
Picco and Pujol,H3 we present a stud y of the interfa- 
cial tension, defined in Eq. ( 3.19 ). High-precision 

data, calculated using L — 2 x 10 8 for 8 < M < 24 on 



the NL very close to the Nishimori point, are shown in 
Fig.|l6|; statistical errors are smaller than symbol sizes. 
As with M/£, one expects, in the critical region and at 
sufficiently large M, to collapse data for M/£ MM onto a 
single curve by plotting it as a function of the scaling vari- 
able (p— p c )M Y l v . Such a collapse is illustrated in FigJlTj, 
using p c = 0.1093 and v = 1.50. Deviations from collapse 
are evident at smaller values of M, appearing as vertical 
offsets of the corresponding lines in Fig.|l7]. Corrections 
to scaling of this type are expected, and arise from scal- 
ing variables which are irrelevant in the RG sense at the 
critical point: in general, we have 



M 2 A/ = a + b{p - p c )M 1/u + cM~ x + 



(4.2) 



where x is the exponent associated with the leading irrel- 
evant scaling variable, a is a universal scaling amplitude, 
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FIG. 14. Data collapse along the NL, using v = 1.50 and 
p c = 0.1093. 
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FIG. 16. Variation of M/^ on the NL close to the Nishi- 
mori point. 



15.0 



10.0 



5.0 



0.0 



•o 




2.33 


□•o 


o M=8 
• M=16 

□ M=32 




\ 


a M=64 


2.23 

=£. 
=L 


\ 




JJLP 






2.13 



-0.20 



-0.10 



0.00 0.10 
(p-pc)M 1/v 



0.20 



FIG. 15. Data collapse along the NL, using v = 1.40 and 
p = 0.1093. 
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FIG. 17. Scaling of M/f w as a function of {p-p c )M 1/u , 
using v — 1.50 and p c = 0.1093. 



and b and c are constants. Such corrections occur at 
the pure Ising transitions and have also been studied in 
the U(l) n etw ork model.a In view of the way that they 
enter Eq. (4.2), it is appropriate to concentrate on the 
M -dependence of the gradients of lines in Fig. 16 when 
determining v. These gradients are shown as a function 
M using a double logarithmic scale in Fig. 18, from which 
we derive our most precise estimate of v, v = 1.50 ±0.03. 

The scaling of M/£ close to the critical point can be 
analysed in just the same way, yielding the same result 
for v. This scaling collapse is depicted in Fig.[l9|. 

We conclude our analysis of critical behaviour on the 
Nishimori line with the results: p c = 0.1093 ± 0.0002 
and v = 1.50 ± 0.03. Our value for p c is consistent with 
the result p c = (U094 ± 0.0002, obtained by Honecker, 
Picco and Pujol,EZl who carried out a detailed study of 



the interfacial tension and correlation functions, using 
the Ising model transfer matrix in a spin basis, which 
restricted system widths to M < 12. Our value for p c is 
also in agreement with some earlier, less precise values, 
including p c = 0.111±0.002, in Ref. ||andp c = 0.1095± 
0.0005 in Ref. |2(| both found using a transfer matrix 
approach with up to 14 spins. It is also marginally in 
agreement with p c = 0.104 from Ref. 29 obtained as the 
critical disorder strength around T = 0. It is in marginal 
disagreement with the result from series expansions ,E£I 
p c = 0.114 ± 0.003. More strikingly, however, our value 
for v is ia disagreement with previous estimates, which 
lie closeES to the percolation value, v — 4/3, including 
most recently v = 1.33 ±0.03 in Ref. ^7|. We believe that 
the larger system sizes accessible in our work, and the 
allowance we have made for irrelevant scaling variables 
at the critical point, together account for the discrepancy, 
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FIG. 18. Scaling of the gradient S of lines in Figjl6| as a 
function of M. The best fit inverse slope is v = 1.50 (solid 
line). Lines corresponding to v = 1.33, 1.45 and 1.55 are also 
shown. 
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FIG. 19. Scaling of M/£ on the NL close the the Nishimori 
point, using v = 1.50 and p c ~ 0.1093 



and that the data shown in Figs. [L5| and |18] exclude this 
smaller value of v. 



E. Scaling along the phase boundary 

The phase boundary separating the ferromagnet from 
the paramagnet coincidesllj with the second relevant scal- 
ing axis at the Nishimori point, in addition to that de- 
fined by the NL. On the boundary, we expect scaling flow 
from C towards the pure critical point for p < p c , and 
from C towards the zero-temperature critical point for 
T < Tn- We analyse such flow in this subsection. 

Qualitative evidence in support of these established 
ideas is presented in Figj2C|, which shows the variation 
of M/£ with position, parameterised by T, on the phase 
boundary, and with M. For p < p c , the coordinates 
of points on the phase boundary are taken from Table 
H while for T < Tn we assume the phase boundary 
to be vertical in the p — T plane and set p — p Cl us- 
ing our estimate for the value of p c . At temperatures 
T > Tn — 0.9533, M/£ decreases with increasing M, 
approaching zero which is the value taken by this scaling 
amplitude in the pure Ising model at T — Tq ~ 2.269; 
fluctuations visible in Fig.|2C| for data at temperatures 
T > 1.5 arise from errors in determining the position 
of the phase boundary. At the Nishimori point itself, 
curves of M/£ for different M cross, with a limiting value 
for M — > oo, as already determined in our study of be- 
haviour on the NL. For T < Tjy, values of M/£ increase 
both with decreasing T and with increasing M, as ex- 
pected if flow is towards lower temperatures. 

Scaling flow along the phase boundary close to the 
Nishimori point is characterised by a critical exponent 
vt, which in principle can be determined using an ap- 
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FIG. 20. Variation of M/£ with position, parameterised by 
T, on the phase boundary. 



proach similar to that taken for v. In practice, there are 
extra difficulties. First, in contrast to the NL, the form of 
the phase boundary is not known exactly; we choose the 
simpler regime, T < Tn, and set p to our estimate for p c , 
as above. Second, it happens that vj- > v, so that flow 
away from the multicritical point is faster in the direc- 
tion of the NL than along the phase boundary. Because 
of this, the range for T over which useful data can be 
collected is limited on both sides. The distance, Tn — T, 
from the Nishimori point should not be too large, or data 
will lie outside the critical region. It should not be too 
small, either, because close to C errors in our value for p c 
will be dominant. Having limited the range for T — Tn 
in this way, the variation in M/£ is also restricted. It is 
therefore particularly important that statistical errors are 
small, and so we study samples of length L — 2 x 10 8 with 
8 < M < 32. The scaled data are presented in Fig.pl]: 
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as with the analysis p rese nted in Fig. fl7| and Fig. |T^, and 
as expected from Eq. 4^2 , the value of vt is determined 
mainly from the gradients of curves for each M. We con- 
clude that vt — 4.0 ± 0.5. While this confidence margin 
is wide, it is encouraging that on extrapolating the data 
in Fig.|Tj to T = Tn we obtain at the Nishimori point 
M/£ = 1.58 ± 0.01 for M — > oo, in perfect agreement 
with the value found independently from data collapse 
on the NL. 
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FIG. 21. Scaling of M/£ on the phase boundary below the 
Nishimori point, using vt = 4.0 and TV = 0.9533 



F. Behaviour at strong disorder 

In three or more dimensions, the random bond Ising 
model has a spin-glass phase at low temperature and 
strong disorder .EJ It is known that spin-glass order does 
not occur in .the two-dimensional RBIM, except at zero 
temperaturej£3 but it is of interest to examine behaviour 
at strong disorder using the methods we have developed. 

Finite-size effects in the RBIM are large at strong dis- 
order and low temperature, as remarked in connection 
with Fig.^, and as is clear from Fig.^||, which shows the 
variation of M/£ with p and M at a fixed temperature, 
T = 0.5, below the Nishimori point. Despite these finite- 
size effects, it is straightforward to identify the position 
of the phase boundary from Fig.^2[ Moreover, the size- 
dependence of M/£ in the paramagnetic phase at T = 0.5 
and higher temperatures is consistent with a finite limit- 
ing value for £ as M — * oo, as required from the fact that 
the RBIM does not have a metallic phase.Q 

For a quantitative analysis of behaviour in this region, 
we focus on the line p — 0.5 which, by symmetry argu- 
ments, is an exact scaling axis. Scaling flow is from the 
zero-temperature fixed point at p = 0.5 towards infinite 
temperature, and one can collapse data on this line to 
extract the limiting behaviour of £ for M — » oo. This 
extrapolated localisation length, £buik, is expected to be 




FIG. 22. Variation of M/£ with p, crossing the phase 
boundary at T = 0.5. 




FIG. 23. Variation of £buik with T on the line p — 0.5: 
ln(£buik) as a function of 1/T. Dashed lines represent 
£buik oc exp(-2/T). 



finite for T > 0. Its temperature dependence for T > 0.4 
(obtained using 8 < M < 64 and L = 10 6 ) is shown in 
Fig. where we compare our, resu lts with the behaviour 
6mik oc exp(-2/T), suggestedBEl for the ±J RBIM. In 
Fig. ^ we compare our same results with the power-law 
divergence, £buik oc T~ v , expected in a RBIM with a 
distribution of bond strengths continuous at J = 0, for 
which exponent values in the range-L' = 3.4 to v = 4.2 
have been reported previouslyoHH Our data in the 
temperature range accessible do not provide firm grounds 
to prefer one form for the temperature dependence over 
the other. 
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FIG. 24. Variation of £buik with T on the line p — 0.5: 
ln(£buik) as a function of ln(T). The dashed line represents 
£buik oc T~ v with v = 4. 



G. Spin-spin correlations 



As a demonstration of the effectiveness of the method 



set out in Sec. Ill C for obtaining even powers of spin-spin 
correlation functions, we have calculated [(erf (n)crj (n)) 2 ] 
at all separations of spins across the width of a long 
system with M = 40. Data at p = 0.08, obtained by aver- 
aging over 10 4 disorder realisations, are shown in Fig.|25|, 
for a high temperature, T = 1.9, lying in the paramag- 
netic phase, and for a lower temperature, T = 1.3, lying 
in the ferromagnetic phase. It is clear for this second 
case that the value of the square of the magnetisation 
can be obtained from the correlation function at separa- 
tions close to Mj 2. 
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FIG. 25. Variation of the disorder-averaged square 
spin-spin correlation function with distance around a sys- 
tem of circumference M — 40, in the paramagnetic phase 
(T = 1.9) and the ferromagnetic phase (T = 1.3). 



We have also used this approach to calculate 
[(erf (n)crj(n)) 2 ] and [(erf (n)erj(n)) 4 ] on the NL at our es- 
timated position for the Nishimori point. At this point, 
one expects decay of the the disorder-average of the 
fc-th power of the spin-spin correlation function to be 
characterised by an exponent r]k- Following the analy- 
sis described in Ref. ^?], and taking M = 20, L large 
and 10 4 realisations, we obtain r\i = 0.183 ± 0.003 pid 
?74 = 0.253 ± 0.003, in agreement with earlier results! 



V. SUMMARY 

To summarise, we have described in detail a mapping 
between the two-dimensional random-bond Ising model 
and a network model with the symmetries of class D 
localisation problems. Building on Refs. we have 
shown in particular how separate boundary conditions 
arise in the network model for sectors of the Ising model 
transfer matrix with even and odd parity under spin 
reversal, and how statistical-mechanical quantities, in- 
cluding the free energy per site and correlation func- 
tions, may be obtained from calculations using the net- 
work model. Amongst other things, this makes clear the 
sense in which the Ising model correlation length may 
be equated with the network model localisation length. 
From a computational viewpoint, calculations based on 
the network model are much more efficient than their 
equivalent using an Ising model transfer matrix in a spin 
basis. This is illustrated by the fact that such calcula- 
tions have in the past mainly been restricted to systems 
of width M < 14 spins, while we present results in this 
paper for M < 256 spins. Applying these ideas to study 
the Nishimori point for the ±J RBIM, we obtain a value 
for the exponent v which is significantly different from 
previous estimates based on much smaller systems sizes; 
our value excludes the possibility of a simple connec- 
tion between behaviour at this critical—point and classi- 
cal percolation, conjectured previously.E3 Beyond compu- 
tational advantages, the equivalence between the RBIM 
and the network model has theoretical interest. It links 
the transition between paramagnet and ferromagnet to a 
version of the quantum Hall plateau transition, as our re- 
sults illustrate. Moreover, even in quasi-one dimensional 
systems for which there is no sharp Curie transition, a 
topological distinction emerges within the network model 
between two separate localised phases. 
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APPENDIX A: EFFECT OF ROUNDING 
ERRORS ON LYAPUNOV EXPONENTS 

The numerical results presented in this paper were ob- 
tained using a modified version of the standard algo- 
rithm for study ing random matrix products, as we de- 
scribe Sec. [II A. The need for such a modification stems 



from the instability of the standard algorithm to round- 
ing errors if the value of the smallest positive Lyapunov 
exponent, e\, approaches zero. The instability is extreme 
and it is of interest to understand how it arises. In this 
appendix we illustrate its origin by examining a simple 
model problem. 

It is sufficient to consider only products of 2 x 2 ma- 
trices, because the instability involves only the space 
spanned by the vectors associated with the pair of Lya- 
punov exponents smallest in magnitude, denoted by 
x M (i) and x M+1 (L) in Sec. Ill A. We therefore consider 
a product of random matrices, each of the form 



T n = 



cosh 9„ sinh 0„ 



sinh 9 n cosh 



(Al) 



and drawn independently from a distribution which has 
(9 n ) = in order that the Lyapunov exponents of the 
matrix product are zero. To model the operation of 
the standard algorithm, we consider evolution of a two - 
com ponent vector v„ under an analogue of Eqs.(3.2)- 
0: 



'n+l 



T„v r , 



(A2) 



In the absence of rounding errors, v„ converges with in- 
creasing n to one of the eigenvectors of T n , and so it is 
natural to expand v„ in this basis, writing 



- 1 
V2 { 1 




1 



-1 



(A3) 



In this notation, Eq.(|A2j) may be written tan0„ + i = 
exp(— 29 n ) tan0„, and has fixed points <p n — rmr/2 
with m integer. We concentrate on the vicinity of one 
of these, considering the range < <f> n <C 1. Then 
fan+i ~ exp(— 29 n )cj) n . We take the effect of rounding 
errors into account by substituting for this the evolution 
equation 



t> n+ i = cxp(-26>„) 



(A4) 



where r\ n is random with (r) n ) = and (r] n r] m ) = 5 mn a 2 . 

A simple treatment of the stochastic process defined 
in this way is sufficient for our purposes. To find ap- 
proximately the limiting distribution P{4> n ) at large n, 



we divide the range under consideration for <p n into the 
regimes < <f> n < a and a < <j) n . In the former the noise 
dominates, generating an approximately uniform distri- 
bution for (j) n . We take 



P{fai) = Cl 



(A5) 



where C\ is a constant. In the latter regime we neglect 
the noise and use in place of <fi n the variable y n = log(0 n ), 
taking its evolution to be 

Vn+l = Vn - 2#n ■ (A6) 

Since we have chosen (9 n ) = 0, this generates a uniform 
distribution for y n in the range a < y n < Y , where the 
upper limit Y ~ log(7r/4) represents the point at which 
the linearisation of tan(</> n ) fails, and also the boundary 
separating the vicinities of the fixed points of Eq.( A2) at 
4> n = and at <j) n = tt/2. On transforming back to <p n 
we obtain within our approximations 



P{4>) 




for < cf) < a 
for a < < 7r/4 



(A7) 



where C2 = C\o for continuity. C\ is determined by the 
normalisation condition 

c7r/4 1 

P(<t> n )dcf> n = - (A8) 



since we may take the full range for fa 
tt/2. We find for a < 1 

d ~ [2(j\n{TT / Aa)]- 1 . 



to be < 



< 



(A9) 



Now consider the effect that noise-induced departures 
of 4> n from the fixed point at <p n = have on the estimate 
of the Lyapunov exponent, e. Using e = (ln|T„v n |), we 
have 

e = i(ln [exp(26»„) cos 2 4> n + exp(-26»„) sin 2 cf> n ] ) . 

(A10) 

Taking for simplicity 9 n and cj> n small, we find 

e ^4(0 2 )(O- (AH) 

In the absence of noise, <p n = and hence e = 0. With 
noise present we must evaluate 

ctt/4 

P{<t> n )<t>ld4> n . (A12) 
1 we find, for a <C 1, 

(A13) 







Using our approximate form for P( 
(fa\) oc I ln(cr)| _1 and hence 

e cx I ln(cr)| 1 



Thus small rounding errors may be responsible for a large 
error in the value obtained for the Lyapunov exponent. 
In the language of this appendix, the modified algorithm 
described in Sec. [II A uses the known symmetry of the 
tranfer matrix to fix <f> n = 0. 
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